% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
%   This scripts is for using the ET ratio to multiply CLM Output model data.
%   It should be run in linux version matlab because of the "/"
%   Dayon 2021/01/21 
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 

clear all
clc
% Read the ET_ratio from nc file
inDir = '/home/dayon/Downloads/GLEAM/ET_ratio';
inFile = 'Global_ETactual_ratio_GLEAM_yr_1986_1999_mo_01_12.nc';
inDirFile = [inDir, '/', inFile];
ET_ratio = ncread(inDirFile, 'ET_ratio');

modDir = '/home/dayon/Downloads/GLEAM/Gleam_Calibration';
% modDir  = 'D:/GLEAM_Monthly/testnc'; 
modFilestruct = dir([modDir, '/','*.nc']);
nFile = length(modFilestruct);
for i = 1:nFile
    modFilename = modFilestruct(i).name
    modDirFile = [modDir, '/', modFilename];
    ET_mod_qsoil  = ncread(modDirFile, 'QSOIL');
    ET_mod_qvege  = ncread(modDirFile, 'QVEGE');
    ET_mod_qvegt  = ncread(modDirFile, 'QVEGT');
    pfts1d_lon  = ncread(modDirFile, 'pfts1d_lon');
    pfts1d_lat  = ncread(modDirFile, 'pfts1d_lat');
    npft  = length(pfts1d_lon); 
    iMonthStr = modFilename(end-13:end-12);
    iMonth    = int32(str2num(iMonthStr));

    parfor j = 1:npft
        indLon = pfts1d_lon(j);
        iLon   = (indLon - 0.5) + 1;
        indLat = pfts1d_lat(j);
        iLat   = (indLat + 89.5) + 1;
        ET_scale = ET_ratio(iLon, iLat, iMonth);
        ET_calib_qsoil(j, :) = ET_mod_qsoil(j, :) .* ET_scale;
        ET_calib_qvege(j, :) = ET_mod_qvege(j, :) .* ET_scale;
        ET_calib_qvegt(j, :) = ET_mod_qvegt(j, :) .* ET_scale;
    end

    % change nc variable value
    ncid = netcdf.open(modDirFile, 'WRITE');
    qsoilid = netcdf.inqVarID(ncid, 'QSOIL');
    qvegeid = netcdf.inqVarID(ncid, 'QVEGE');
    qvegtid = netcdf.inqVarID(ncid, 'QVEGT');
    netcdf.putVar(ncid, qsoilid, ET_calib_qsoil);
    netcdf.putVar(ncid, qvegeid, ET_calib_qvege);
    netcdf.putVar(ncid, qvegtid, ET_calib_qvegt);
    netcdf.close(ncid);
end